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Abstract 

A solution to the ultra-relativistic strong explosion problem with a 
non-power law density gradient is delineated. We consider a blast wave 
expanding into a density profile falling off as a steep radial power-law 
with small, spherically symmetric, and log-periodic density pertur- 
bations. We find discretely self-similar solutions to the perturbation 
equations and compare them to numerical simulations. These results 
are then generalized to encompass small spherically symmetric per- 
turbations with arbitrary profiles. 
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1 Introduction 



The profusion of explosions occurring in the observable universe has led to a 
pointed interest in the dynamics of blast waves. The quantitative treatment 
of strong explosions began with the self-similar solutions found by Sedov, 
Von-Neumann and Taylor [1] [2] [3] for the flow behind spherical Newtonian 
shocks propagating into a cold gas with a power-law density profile, p oc r~ k . 
This solution is valid for moderately steep decay exponents, k < 3. Subse- 
quently, corresponding solutions were found in the ultra-relativistic regime by 
Blandford & McKee [I] for k < 4. In these solutions the Lorentz factor of the 
shock T scales as T 2 oc t~ m , and m is fixed by energy conservation arguments. 
For k > 4 this procedure fails due to the energy in the Blandford-McKee solu- 
tion diverging, and a different argument must be used if a self-similar solution 
is to be found. It turns out that there exists a k g > 4 such that for k > k g just 
such an argument exists [5]. The solutions in this regime are called type-II 
solutions [6J, and what sets them apart from the solutions with k < 4, known 
as type-I solutions, is the rapid acceleration of the shock and the fluid be- 
hind it that causes the formation of a sonic point between the shock and 
the center of the explosion. This point (actually a spherical surface) marks 
the boundary of an inner region that becomes causally disconnected from 
the shock. This allows a self-similar flow behind the shock to coexist with a 
non self-similar flow further away from it, thus maintaining a finite amount 
of total energy. For the ultra-relativistic case, k g = 5 - v/374 « 4.13. The 
sonic point appears as a singularity in the hydrodynamic equations, and the 
requirement of regularity in traversing this singularity supplies the necessary 
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condition to fix the value of m. 

In this paper we focus on Type-II ultra-relativistic solutions with k > k g 
[S], and use these as the basis for a perturbative analysis where we disturb 
the external density profile. The basic method for doing so was developed 
in a previous paper [7j for the Newtonian case, and we adapt it here for the 
ultra-relativistic case. We introduce a special family of perturbations to the 
ambient medium surrounding the explosion, so that we are able to reduce 
the hydrodynamic equations to a set of ordinary differential equations, and 
then find the flow behind the shock in the presence of perturbations. In 
section [2] we briefly describe the unperturbed solutions, while in section [3] 
we shed light on some general properties of these solutions. In section H] 
we write down the equations for the perturbations and solve them, and in 
section [5] we generalize the results of the previous section by using a spectral 
decomposition of an arbitrary perturbation profile. Finally in section [H] we 
summarize and discuss our results. 

2 The Unperturbed Solution 

We begin with a quick review of the self-similar solutions that will later serve 
as the basis for perturbation. These describe the flow after the discharge of 
a large amount of energy in a small volume surrounded by a spherically 
symmetric distribution of stationary cold gas. The density of this medium 
follows a radial power law, 
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A more detailed treatment is given in Best & Sari [5j. We use units such 
that the speed of light is unity throughout the paper. 



2.1 The hydrodynamic equations 

The hydrodynamic equations for a relativistic ideal gas are derived from the 
local conservation of energy, momentum and number of particles. In spherical 
symmetry they take the form 
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where the velocity /3 and the particle number density n' are measured in the 
frame of the unshocked gas, and the pressure p and internal energy density 
e are measured in the fluid rest frame. The particle density in the fluid 
rest frame is related by n! = to that in the explosion frame, with 7 = 
(1— f3 2 )~ l l 2 . We will treat these equations in the ultra-relativistic limit where 
T >> 1, and use the relativistic equation of state: 
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We take the Lorentz factor of the shock to scale as 



T 2 oc r 



(4) 
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so the shock radius is, to first order in 1/r 2 , 



R(t) = t 



1 - 



2(m + l)r 2 _ 



(5) 



where m is as yet undetermined. We now make a change of variables to a 
similarity coordinate that follows the scale height in the problem, R/T 2 . The 
coordinate x is defined by 



X = [l + 2(m + l)r 2 ](l-r/t). 



(6) 



We define the self-similar functions /, g and h as 



P = ^wiT 2 f{x) 
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n = 2niT 2 h(x) 



(7) 



where we can take wi, the enthalpy before the shock, equal to p\ because the 
unshocked gas is assumed to be cold. 

Substituting equations (J6j) and (JTj) into equation (j2J) we get the hydrody- 
namic equations in their self-similar form: 
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The boundary conditions for these equations are derived from the Rankine- 
Hugoniot conditions at the shock which can be written as 



2^2 

P2 = -r wi 



n' 2 = 2r 2 ni 

il = ^r 2 , (9) 



where subscript l's denote quantities just ahead the shock and 2's quantities 
immediately behind the shock. In terms of the self-similar quantities, these 
boundary conditions become 

f(x = i) = g(x = i) = h( x = i) = i. (io) 

The self-similar form (jSJ) of the hydrodynamic equations exposes the sin- 
gularity at the sonic point, g(Xs)Xs = 2(2 — v3)> where the denominators 
vanish. A valid solution should be regular everywhere, including the sonic 
point, and so we can require the numerators as well as the denominators to 
vanish there. This condition yields the proper value of m, 



m 



k(3 — 2y/3) — 4(5 - 3\/3) (IT 



The equations (jSJ) may be solved implicitly in terms of the variable 



3X + W3-2y^fc-4 
X 10^3-2^3*; -3 
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to give 

log(g) = (k-Q)(3-2V3)log(x) 

log(f) = (k-Q)(A-2V3)log(x) 

(k - 6) [33 - 20^ + (4^ - 6)k)log{x) + (k + 2^ - A)log(2 - g X ) 

l-5V3 + V3k 

(13) 



log(h) 



These solutions now form the basis of the perturbative analysis to be dis- 
cussed in section HI 



3 Type-II solutions as simple waves 

The hallmark of Type-II flows is the existence of a sonic point. We now 
digress to explore a special property that arises from the combination of a 
sonic point with ultra relativistic self similar flows. We begin by considering 
the Riemann invariants (RI) of the flow [8j. Simply put, the RI are quan- 
tities conserved along characteristics of the flow. Considering only the C± 
characteristics here, this condition may be written as an advection equation 
for the respective invariants J±: 



d u±c d 

777 + 



dt 1 ± uc dr 



J± = D ± J ± = 0. (14) 
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where c is the local speed of sound and D is defined as the advective deriva- 
tive. For the case of plane waves, the RI are given by [9] 



jplanar = I dp ± ( l+£ ^ 



c{e + p) 2 \l-f3J' 
Which evaluates under the relativistic equation of state ([3j) to 



Ji anar = ^ln(p)± l -ln^). (16) 



However, it may be seen by substitution in equation (j2J) that for spherical 
geometry these quantities satisfy 



D jplanar = /^X 

± r 1 + uc y ' 



The right hand side of (117j) must be incorporated into the advective term 
on the left for us to find the spherical J± that satisfy (fT%l) . This is in general 
fruitless, but we may here take into account the nature of relativistic self- 
similar flows, and find a solution valid under the ultra-relativistic approxi- 
mation prevalent throughout this work. As we saw earlier, the solutions at 
hand are concentrated on a thin shell with thickness of order R/T 2 . That 
means that in the limit of high Lorentz factors we may safely take r pa R pa t, 
and substitute that into (fTTj) . Further taking /3 pa 1 and c = l/v3, we arrive 
at 

D ± ^ anar + —l—lnit^j = 0. (18) 
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so we may write 

J± = ^ p ±V3/2 t ±mV3+i)_ (19) 

Plugging in the solutions (fT3l) into ffT9l) reveals that J + is a constant 
(both along C + characteristics and across them). This is somewhat surprising 
because there is no reason for either of the RI to be constant in a general 
flow. We can understand how it comes about by considering the special 
causal structure created by the presence of a sonic point. The sonic point 
acts as a stationary point for outgoing C + characteristics, not unlike the event 
horizon for outgoing light rays in a gravitational black hole. Thus we can 
see , as shown in figure [TJ that all C + characteristics asymptotically approach 
the sonic point when traced back to early times. This means that the C + 
characteristics that eventually fan out to cover all the space below the shock 
carry with them a single value of J + , the one at the sonic point. 

Flows with a constant value of one RI are called simple waves. The 
crucial consequence of this property is that sound waves traveling along the 
direction opposite to that associated with the constant RI (ingoing waves in 
our case) will do so without being scattered, even though they are traversing 
a non uniform and non steady background. This can be understood by 
considering the role of RI in wave scattering. If we consider (without loss of 
generality) left going sound waves in a general flow, these waves will create 
perturbations to the fluid velocity and to the speed of sound. These in turn 
will cause the right going characteristics to be slightly deflected, and thus the 
value carried by J + to some point P, as depicted in figure [21 will deviate from 
the same value without the left going waves, giving rise to scattered right 
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going waves. However, if J + is constant, and the hydrodynamic variables are 
definite functions of the RI, this has no effect. The same value of J + will 
arrive at P, and no change in the hydrodynamics will be felt there. 

The latter condition is satisfied in our case by virtue of the relativistic 
equation of state. If a general relation exists between p, e and p the integral 
in (|15|) will be path dependent and thus hydrodynamic variables at P will 
depend both on the value of J+ and on the history of the C + characteristic 
that arrived there, nullifying our previous result. We now come to the con- 
clusion of this section, where we see that several effects conspired, in the case 
discussed here, to make our flow a simple wave and thus greatly simplify the 
perturbation analysis. 

4 Perturbations 

Applying an arbitrary perturbation to the ambient density profile ([T]) will 
in general lead to a non self-similar flow that will be difficult to solve for 
without a full numerical simulation. We can however carefully choose a 
special perturbation that results in tractable equations for the perturbed 
flow. We consider a log-periodic perturbation, 



The merit of this choice is that the perturbation wavelength scales like the 
radius, introducing no new scale into the problem. The parameter a is a 
small dimensionless amplitude, and (3 is the logarithmic frequency. r has 




(20) 
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dimensions of length and only serves to determine the phase of the pertur- 
bation. While the physical perturbation is only the real part of 5p in ( 1201 ). 
writing it as a complex power law facilitates finding a self-similar solution. 
We rewrite the perturbed hydrodynamic functions and shock radius as 



p =l Wl r 2 [f( x ) + b(t)5f(x)} 



l 2 = l^[9(x) + b(t)5g(x)} 



n> = 2n 1 T 2 [h(x) + b(t)5h(x)] (21) 



and 



R + 5R = t 
where we take b(t) to be 



1 - 



2(m+ l)r 2 



'.22) 



b(t) = a -R q . (23) 

We now have two unknown parameters: q, that describes the frequency of 
the perturbations behind the shock, and d which describes their complex 
amplitude. Of these, q may be readily seen from the boundary conditions 
to satisfy q = i/3. As can be expected in a linear perturbation analysis 
this means that the perturbations behind the shock oscillate at the same 
frequency as the disturbance before the shock. The value of d, which connects 
the strength of the perturbations ahead of and behind the shock, is less 
obvious and must be solved for along with the perturbations. We can now 
substitute ( 12TI) into (j2J) and linearize with respect to b(t) to obtain equations 
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for Sf, Sg and Sh. These equations may be written as 



MY' + LY = 0, 



(24) 



where Y(x) = (Sf, Sg, Sh)" J 
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It is notable that the equations for the pressure and velocity are decoupled 
from the density, a consequence of the relativistic limit where the rest mass is 
a negligible part of the energy. These equations must be solved in conjunction 
with the proper boundary conditions. The conditions at the shock are derived 
from equation (JH]) by comparing the perturbed functions at the location of 
the unperturbed shock with the unperturbed functions at the same point, to 
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first order in b. This yields 



Sf(x = i) 



q + m + 1 



+ d + f(l) 



Sg(x = 1) 



m + 1 
q + m + 1 




M(X = 1) 



m + 1 
q + m + 1 



(27) 



m + 1 



One additional condition is needed to close this system and determine the 
value of d, and that is the boundary condition at the sonic point. The 
condition of regularity at the sonic point is again satisfied if the numerators 
in equation (j2~"Ij) vanish at the sonic point. This yields the condition 



which together with equations (124p and (1271) enables us to find both d and 
the functions 5f, 5g and Sh for any k and (3. The straightforward solution 
would be to start with (12*71) at the shock and use a shooting method to find 
the value of d that satisfies equation ( 1281) at the sonic point. 

We can, though, use the insight provided by the discussion in section [S] to 
simplify the problem. The shock acts as a source of ingoing sound waves that 
propagate as the perturbations to u and p that we are seeking. As we have 
shown in section [5] the unperturbed solutions is a simple wave, and so these 
ingoing waves travel without scattering towards decreasing radii. Beyond 
the sonic point all characteristics point at decreasing values of x, and thus 
no reflections may travel back out. This is true wherever the self-similar 
solution is valid, which is our scope of interest. This region covers a fraction 



2 Sg(xs) 




f(Xs) 



Vs g(xs) 
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of the volume inside the shock that asymptotically approaches unity. We 
thus see that to solve the perturbations equations we need only take into 
account ingoing waves. The form of J + may be differentiated to reveal the 
connection between pressure and velocity perturbations: 

P \/3 7 2 

Equation (j28p therefore holds everywhere and not only at the sonic point. 
This relation enables us to find 5g, Sf and d analytically: 



Sg = (2V3 + 0) 2-2 m -t + ,:/3V 1 +^ x -U- 6j+m _, fi 
m + 1 ^4 - 8g X + g 2 x 2 

<~(l + V>(l + J^). (30) 



We are left only with the equation for 5h that may be solved using numerical 
methods. The value of d reveals that at large frequencies the amplitude of 
perturbations behind the shock is inversely proportional to the frequency of 
the density perturbations (3, and the phase is a quarter wave behind them. 
Representative solutions are shown in figure E] along with corresponding re- 
sults of a full numerical simulation performed with a second order Godunov 
type scheme. 

It bears mention that while /, g and h are functions of x only, the solution 
we find is not strictly self-similar. It is rather composed of one truly self- 
similar part, the unperturbed solution, and another part which is discretely 
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self-similar, the perturbation. While the perturbation appears to be in itself 
continuously self-similar, this is only because we used complex notation to 
write it. Once we take the real part of the complex exponential in (123]) the 
oscillatory nature of the perturbations is revealed and we get only discrete 
self-similarity, with a period of 



5R a* 

— = e--l, (31) 



very much like the Newtonian case treated in Oren & Sari [7j. 

5 Arbitrary Perturbations 

The method outlined so far is limited by the special choice of density per- 
turbations. However, since these special perturbations form a complete set 
we can use them to decompose any given perturbation. The linearity of our 
scheme ensures that summing the solutions with the appropriate weights will 
give us a solution for the compound perturbation, in a similar fashion to the 
well known Fourier method for solving differential equations. The resulting 
solution will no longer be self-similar, due to the different time dependence 
of each component. We demonstrate the validity of this procedure by consid- 
ering a scenario commonly encountered in astrophysical context, an abrupt 
density jump: 

p(r) = Kr~ k [l + s9(r - rj)] (32) 
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where 6(x) is the Heaviside step function, is the location of the discontinu- 
ity and s is a small amplitude. In figure H] we compare the solution obtained 
by a fourier decomposition of the step function to a full numerical simula- 
tion with the corresponding initial conditions, at a specific point in time. 
The discontinuity depicted there is a reverse shock propagating back into the 
shocked material and slowing it down. 

6 Discussion 

We have laid out a method of solving the hydrodynamic equations for a 
strong explosion in the presence of spherically symmetric perturbations to 
the ambient density. At first we study a special group of perturbations with 
log-sinusoidal radial dependence, and discover an analytic solution. The 
requirement of spherical symmetry is not easily relaxed, because relativistic 
effects make it difficult to find self-similar solutions with a non trivial angular 
dependence. Another limitation is the linearity of the perturbation analysis 
that limits the validity of the solutions to small amplitudes. The smallness 
required may be seen in figure where numerical solutions with different 
amplitudes are superimposed against the analytical solution. It can be seen 
that our theory gives reasonably accurate results for amplitudes up to about 
0.1. The nonlinearity of waves with higher amplitudes is expressed through 
shock formation before their crests, as can be seen in the red line in figure El 
On the other hand we take advantage of linearity to generalize our results 
using a Fourier-like method, decomposing an arbitrary perturbation to simple 
modes for which we can solve the equations analytically. In this way we 
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can treat interesting scenarios like a sudden rise or drop in density, which 
might e.g. be encountered in a stellar wind due to interaction with the 
interstellar medium. Another possible application is the emergence of a shock 
from the edge of a star |10j . where the drop in density accelerates shocks to 
relativistic velocities. The effect of perturbations in the star's envelope can 
be treated with the method presented here, requiring only an adaptation of 
the unperturbed solution. 
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Figure 1: The causal structure of Type-II solutions. At late times all of space 
is covered by C + characteristics that originate at the sonic point, making J + 
constant. 




Figure 2: The two sets of characteristics and their interaction. Solid and 
dashed lines represent perutbed and unperturbed characteristic curves, re- 
spectively. In this example perturbations to J_ create a deflection of the C + 
characteristics, thus causing a different characteristic to reach the point P. 
This however will not make a difference at P if J + is constant. 
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Figure 3: The solution to the perturbation equations against the result of a 
full non-linear numerical simulation. The solid black line is the theoretical 
value calculated from equation ( l30l) . and the dotted lines in color represent 
the numerical solution, converted to self similar form using equation (l2"Tl) . 
The four lines represent four consecutive periods of the external perturbation, 
demonstrating the dicretely self-similar nature of the solutions. 
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Figure 4: The relative Lorentz factor perturbation following a sudden in- 
crease in the ambient density, with k = 5, s = 0.05, rj = 10 and R/Tj = 1.75. 
The numerical solution is smoothed to reduce numerical noise. 
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Figure 5: The pressure perturbation for k = 5, (3 = 20 and different values of 
a, normalized by o to facilitate comparison. The solid black line represents 
the analytic solution. 
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